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Abstract 

A  computer  program  for  simulating  coupled 
phenomena  of  propulsion-generated,  chemically 
reacting,  two-phase  jet  flow  fields  is  described  in 
detail.  A  sample  calculation  is  performed  using  typ¬ 
ical  solid-propellant  rocket  motor  operating  condi¬ 
tions  as  input.  The  model  contains  approximations 
for  turbulence,  particle  drag,  chemisty,  and  particle 
supercooling.  The  particle  phase  change  is  simu¬ 
lated  kinetically  and  includes  two  solid  Al203  crys¬ 
talline  structures:  meta  stable  gamma  and  stable 
alpha  forms. 

Introduction 

Accurate  modeling  of  multiphase  alumina- 
loaded  rocket  exhaust  flow  fields  is  essential  for 
assessments  of  the  missile  base  heat  transfer  to 
ensure  the  motor’s  health  and  survival.  Two-phase 
gas/particle  flow-field  simulation  is  also  an  impor¬ 
tant  element  of  efforts  to  predict  solid-propellant 
rocket  motor  effluents  for  environmental  studies 
including  space  debris.  In  order  to  assess  the  nom¬ 
inal  performance  of  solid-propellant  propulsion  sys¬ 
tems,  two-phase  simulations  are  essential.  The 
condensed  particulates  stabilize  the  combustion 
chamber  process  and  increase  the  combustion 
chamber  temperatures;  however,  the  drag  intro¬ 
duced  by  the  particulates  in  the  expanding  nozzle 


represents  a  two-phase  flow  performance  loss  and 
a  decrement  to  the  thrust  chamber  performance. 
Thus,  the  presence  of  particulates  in  thrust  cham¬ 
bers  must  be  carefully  evaluated  in  performance 
assessment  methodologies.  In  reality,  this  problem 
presents  a  scenario  in  which  a  continuous  range  of 
solid,  liquid,  and  multiple  phase,  i.e.  liquid,  as  well 
as  multiple-phase  solid  crystalline  particulates,  are 
spewed  out  within  a  hot  bath  gaseous  mixture  from 
the  combustion  chamber.  Expansion  of  the 
exhaust  gas  acts  as  a  spectator,  tending  the  small 
particles  to  a  higher  velocity,  and  the  larger  parti¬ 
cles  to  a  higher  internal  energy  (temperature).  It  is 
not  surprising,  then,  that  the  implementation  of 
known  particle  phase  change  kinetics  presupposes 
an  intrinsic  coupling  with  the  most  up-to-date  mul¬ 
tiphase  flow  models. 

Many  standard  models  approach  the  problem 
using  an  equilibrium  process  to  describe  the  phase 
change  of  thrust  chamber-generated  aluminum 
oxide  (Al203)  particulates  from  liquid  to  solid  as  the 
flow  expands  and  cools.  In  this  approach,  the  liquid 
particle  temperature  is  tracked  until  the  melting 
temperature  is  reached.  At  this  point,  an  equilibrium 
phase  change  process  is  initiated  where  the  heat  of 
fusion  is  released  from  the  particulate  to  the  lower 
temperature  gas  stream.  The  particulate  tempera¬ 
ture  is  maintained  at  the  melting  temperature  until 
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the  heat  of  melting  is  completely  absorbed  by  the 
gas  stream.  At  this  point  the  particulates  are  in  the 
stable  solid  alpha  phase.  In  an  expanding  flow  envi¬ 
ronment,  the  temperature  of  the  solid  particle  will 
continue  to  decrease  as  the  gas  flow  continues  to 
expand  and  cool.  This  phase  change  model  is  an 
oversimplification  because  it  does  not  account  for 
the  crystallization  kinetics  and  the  metastable 
gamma  crystallization  phase  which  have  been 
experimentally  observed  in  the  laboratory  during 
analyses  of  the  phase  change  process.  We  have 
taken  an  alternate  approach  and  have  solved  the 
equations  of  motion  for  all  phases  of  the  exhaust 
flow,  including  gas,  liquid  particles,  particles  of  both 
liquid  and  solid  phases,  i.e.,  slush  balls,  and  parti¬ 
cles  of  multiple  solid  crystalline  structures.  The 
most  significant  advancement  of  this  work  relative 
to  the  industry  standard  is  that  the  particle  phase 
change  process  is  kinetically  controlled  and 
includes  the  transition  from  metastable  gamma 
crystalline  structure  to  the  stable  alpha  phase.  This 
model  also  retains  the  essential  physical  character 
of  particulate  undercooling  also  observed  in  the 
laboratory. 


Mathematical  Models 

Gasdynamics 

The  system  of  gasdynamic  equations  for  planar 
(v  =  0)  or  axisymmetric  (v  =  1 )  flows  can  be  written 

in  the  following  form: 

3x  3y  y 
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Here  x  and  y  are  the  axes  of  a  Cartesian  system 
of  coordinates,  u  and  v  are  the  velocity  components 
in  these  coordinates,  p  and  P  are  the  density  and 
the  static  pressure,  h0  is  the  specific  total  enthalpy, 
qn  is  an  additional  parameter  which  describes  the 
nonequilibrium  processes  occurring  in  the  gas  mix¬ 
ture,  $  is  a  term  defined  by  the  kinetic  theory  of 
nonequilibrium  processes,  and  F  is  the  term 
responsible  for  viscosity,  heat  conductivity,  diffu¬ 
sion,  and  interaction  with  the  particles. 


Chemical  Reactions 

One  assumes  that  there  are  several  indepen¬ 
dent  chemical  reactions  in  a  gas  mixture.  The  form 
of  the  gas  phase  chemical  reactions  is  as  follows: 

yc+  a_  =  yc  a„ 

Z-i  mn  n  mn  n 

n  n 

where  C^n  and  C^n  are  stoichiometric  coeffi¬ 
cients  specifying  the  quantity  of  An,  the  chemical 
component  taking  part  in  the  forward  and  reverse 
directions  of  the  mth  reaction,  respectively.  Thus,  as 
the  result  of  mth  forward  or  reverse  reaction, 
ACmn  =  C~n-C*n>  molecules  of  An  type  are 
formed  or  destroyed. 


For  chemical  kinetic  equations,  the  mole-mass 
concentrations,  an,  are  chosen  as  qn.  The  terms 
3>n  in  the  right-hand  side  of  Eq.  (5)  have  the  form 


= 


RT 


n  PM, 


£AC mn  >dI<P»-/RT> 
m  L  n' 


mn 


-  K;ll<pn-/RT> 


where  T  is  the  static  temperature,  M£  =  1/£an  is 
the  molecular  weight,  R  is  the  universal  gas  con¬ 
stant,  Pn  =  rnP  is  the  partial  pressure,  rn  =  Mxan  is 

the  mole  fraction,  and  K  +  (T)  and  K  (T)  are  the 

m  m 

rate  constants  for  the  forward  and  reverse  direc¬ 
tions  of  the  mth  reaction,  respectively. 


2 

American  Institute  of  Aeronautics  and  Astronautics 


Table  1  includes  the  chemical  rate  equation  and 
a  set  of  chemical  reactions  and  coefficients  used  to 
determine  the  rate  constants  for  combustion  prod¬ 
ucts  of  the  C-H-O-N-CI  type.  They  are  chosen  in 
accordance  with  Ref.  1 ,  Vol.  1  and  Ref.  2,  Chapter 
17.  The  reverse  direction  rate  constants,  Km(T), 
are  defined  using  the  equilibrium  contants 
Km(T)  =  K*(T)/K~(T)  which  are  taken  from  Ref.  3. 


Laminar  and  Turbulent  Mixing 


For  predictions  of  the  viscous  flow,  the  terms 
Fvis  must  be  added  to  the  right-hand  side  of  Eqs. 
(2)-(5).  In  the  parabolized  Navier-Stokes 
equations4  (where  streamwise  diffusive  terms  are 
eliminated),  these  terms  can  be  written  as 


_vis  _  3U  U  pVis  _  3V 
u  "  3y +Vy’  v  “  ity 

F;is  =  |i+vt!,Fvnis3+v^ 

h  3y  y  n  3y  y 


U  = 


-K 


/  ,3v 

<H  +  m>3y 


(6) 


Here,  p  =  p(T,  q)  and  pt  are  the  laminar  and  tur¬ 
bulent  viscosity  coefficients,  h  =  h(T,q)  is  the  spe¬ 
cific  static  enthalpy,  and  Pr  and  Prt  are  the  laminar 
and  turbulent  Prandtl  numbers.  In  Eq.  (6)  one 
assumes  that  the  laminar  and  turbulent  Schmidt 
numbers  are  equal,  Sc  =  Pr  and  Sct  =  Prt  =  0.7. 
The  term  (1  +  v2/u2)  is  inserted  into  the  equations 
to  take  account  of  the  shear  layer  inclination  with 
respect  to  the  x  axis. 


In  turbulent  flow  analysis,  u,  v,  p,  and  P  are 
treated  as  the  Favre-averaged  values.  A  variant  of 
the  K-e  turbulence  model  suggested  in  Ref.  4  is 


Table  1.  Chemical  Reactions  and  Rate  Constants 


K+  =  AT  n  exp(-E/RT) 


Reaction 

A,  m6/(kmol-sec)2 

n 

E,  cal 

GO  +  0  +  M  =  CO2  +  M 

3.5  •  10° 

0 

2100 

OH  +  H  +  M  =  H20  +  M 

1.2  -  1014 

1 

0 

0  +  N  +  M=  NO  +  M 

3.3-  109 

0 

0 

h  +  h  +  m  =  h2  +  m 

1 .4  -  1014 

1.5 

0 

o  +  o  +  m  =  o2  +  m 

5.5-  1011 

0.87 

0 

N  +  N  +  M  =  N2+  M 

2.7  •  1010 

0.5 

0 

H  +  0  +  M  =  OH  +  M 

3.3  •  1012 

0.5 

0 

Cl  +  Cl  +  M  =  Cl2  +  M 

3.6  •  108 

0 

-1800 

H  +  Cl  +  M  =  HCI  +  M 

1.4- 1016 

2 

0 

Reaction 

A,  m3/kmol-sec 

n 

E,  cal 

OH  +  CO  =  C02  +  H 

2.5-  109 

0 

5100 

H2  +  OH  =  H20  +  H 

1.1  •  1011 

0 

8600 

OH  +  OH  =  H20  +  0 

1.0  -  1010 

0 

1200 

H2  +  0  =  OH  +  H 

1 .3  •  1010 

0 

9860 

02  +  H  =  OH  +  0 

2.2  •  1011 

0 

16,500 

02  +  N2  =  NO  +  NO 

5.2  •  1010 

0 

107,000 

NO  +  N  =  N2  +  0 

3.0-  1010 

0 

200 

NO  +  0  =  02  +  N 

1.1  •  1010 

0 

41,700 

Cl  +  H2  =  HCI  +  H 

8.4-  109 

0 

4260 

Cl  +  H20  =  HCI  +  OH 

3.0-  1010 

0 

17,600 

HCI  +  0  =  OH  +  Cl 

3.6  •  1010 

0 

6000 

H  +  Cl2  =  HCI  +  Cl 

9.0  -  1010 

0 

1200 

used  for  pt  determination.  Equations  for  the  kinetic 
energy  of  turbulence  (k)  and  its  dissipation  (e)  coin¬ 
cide  with  Eq.  (5)  if  F  and  $  are  defined  by 


vis  3K  K  Fvis  =  3E  E 
k  3y  y’  e  3y  y 
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-  e(1  +  C^ 

2 

-e(C2  +  C1 ) 


where 


C1  =  min^1.32,  0.3Mt  +  7.5M2^,  Mt  =  Jup/yP 
C2  =  1.92  -0.0667f 

Here,  y  is  the  specific  heat  ratio  and  f  is  set  to 
zero  in  the  region  extended  axially  from  the  nozzle 
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exit  to  the  cross  section  where  the  shear  (mixing) 
layer  reaches  the  x  axis.  Downstream  of  that  cross 
section  location,  f  is  defined  by 

f  =  Vi  duo 02 
(ue-u0)dx 

where  u0  and  u£  are  the  gas  velocities  at  the  x  axis 
and  the  external  flow  velocity,  respectively.  Turbu¬ 
lent  viscosity  coefficient  is  defined  by 

„  k2[1  -0.222(1  +  C^] 

***  ”  C^P  e  0.78 

where  =  0.09  -  0.04f 

Thermodynamics 

The  specific  total  enthalpy  is  defined  by 
h0  =  (u2  +  v2)/2  +  k  +  h,  h  =  ^anhn 

n 

where  hn  =  hn(T)  is  the  mole  enthalpy  of  the  nth 
chemical  component,  including  its  formation  heat. 
These  enthalpies  are  taken  from  the  reference 
handbook.3 


The  specific  heat  ratio,  y,  for  gaseous  mixtures 
is  understood  to  be  "frozen"  and  defined  by 


v  =  IP  c  =  ah 


.cv  =  cD-i^ 


<V  p  3T|q  =  const’  V 
The  equation  of  state  has  the  form 


However,  the  k-e  turbulence  model4  requires 
that  the  additional  pressure  of  turbulent  pulsation 
be  included  in  the  momentum  Eqs.  (2)  and  (3). 
Because  of  this,  the  value  P  will  be  understood  as 
the  sum  of  both  the  static  and  additional  pressures 
created  by  turbulent  pulsations.  Therefore,  the 
equation  of  state  should  be  rewritten  in  the  form: 


p  =  eBI  +  P,  P  =  |Pk[1  - 0.222(1  +^)] 
Mv  3 


Liquid/Solid  Phase  (Particles) 

The  mathematical  modeling  of  multiphase  flows 
in  the  context  of  multitemperature  and  multivelocity 
simulations  must  invoke  additional  equations 
describing  the  motion  of  AI2O3  particles,  including 
the  particle  phase  thermodynamics  and  kinetics. 
Under  the  assumptions  that  the  particles  are 
spherical,  chemically  inert  to  the  gaseous  phase, 
and  do  not  interact  with  each  other,  these  equa¬ 
tions  can  be  written  as: 

du.  dvs 

=  Cf(u-us),-^  =  Cf(v-vs)  (8) 

"3T  *  cq<T_Ts)  <9> 

where 

Cf  =  Cd  =  Cd(Resj  Mg  Y>  Ts/T) 

8rsPs 

C  _  3NuX,  Nu  _ _ N H! _ 

q  2p°r2’  1 +3.42Nu°Ms/ResPr 

Nu°  =  2  +  0.459Re°'55Pra33,  Ms  =  ^ 

2W<5prq  /  2  2 

Res  =  — Ws  =  J(u  -  uf  +  (v  -  vs)2 

Here  index  's’  denotes  the  number  of  particle 
fractions  (all  particles  from  each  fraction  are  of  the 
same  size),  us  and  vs  are  the  particle  velocity  com¬ 
ponents,  d/dt  s  us3/9x  +  vs9/3y  is  the  total  deriva¬ 
tive  with  respect  to  the  time  (along  the  particle  tra¬ 
jectory),  es  and  Ts  are  the  specific  internal  energy 
and  the  temperature  of  particles,  respectively,  rs  is 
the  particle  radius,  p°  =  3329  kg/m  is  the  intrin¬ 
sic  density  of  alumina  material,  X  and  p  are  the 
heat  conductivity  and  the  viscosity  of  the  surround¬ 
ing  gas,  and  a  =  7y(P/p)  is  the  gas  sound  veloc¬ 
ity.  Unitless  drag  coefficient  of  particles,  Cd,  is 
defined  by  relationships  given  in  Ref.  5. 
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Nonequilibrium  crystallization.*  As  Al203 
particle  temperatures  in  the  nozzle  and  exhaust  jet 
decrease,  they  experience  the  liquid-to-solid  phase 
transition  process  releasing  the  crystallization 
energy.  This  is  a  nonequilibrium  process  which  is 
described  using  the  relative  crystallization  front 
radius,  r* .  A  particle  entirely  in  the  liquid  state  is 
denoted  by  r*  =  1 ,  whereas  r*  =  0  corresponds 
to  the  complete  solid  state.  If  the  particle  is  in  the 
transitional  state,  then  0  <  r*  <  1 ,  and  the  kinetic 
recession  of  the  crystallization  front  is  described  by 
the  equation 

drs  1  8 

r»dT  =  -a*<Tm-TS>  <10> 

where  a*  =  0.64  -10“6  m/sec  •  K1-8  and  Tm  =  2289 
K,  the  equilibrium  melting  temperature.  According 
to  the  adopted  model,  particle  crystallization  begins 
just  below  the  critical  temperature,  Tc  =  0.83Tm 
(1900  K).6,7  The  liquid  phase  is  assumed  to  trans¬ 
fer  initially  to  the  metastable  solid  y  phase  for  all 
conditions.  The  transition  process  from  the  y  phase 
to  the  stable  a  phase  begins  as  soon  as  any  portion 
of  the  y  phase  has  appeared.  To  determine  the  por¬ 
tion  of  solid  substance  in  the  a  phase  (fa),  the  fol¬ 
lowing  equation  is  applied8: 

df 

-gf  =  aaexp(-ba/Ts)  (11) 

where  aa,  ba  are  empirical  coefficients.  According 
to  measurements9  and  recommendations,8  ba  = 
58,368  K.  The  value  of  aa  is  chosen  to  be  equal  to 
1.5  •  1012  1/sec.10 

While  the  solid  phase  portion  of  the  particle  is 
known,  the  total  portions  of  the  y  phase  (My)  and 
the  a-phase  (Ma)  are  defined  by  formulae 

My  =  (1-(r;)3)1-fa).Ma=(1-(r;)3)a 

The  specific  internal  energy,  es,  is  connected 
with  the  particle  temperature,  Ts,  by  the  equation 


‘The  mathematical  description  of  Al203  particle  kinetics  has  been 


es  —  CSTS  +  ^s^rsj 

where  hs  =  0.915  •  1 06 J/kg  is  the  latent  heat  of 
the  Al203  fusion  and  Cs  =  1600  J/kg  •  K  is  the  spe¬ 
cific  heat  capacity  of  particles.  Here,  the  coefficient 
Cs  is  assumed  to  be  equal  for  both  the  liquid  and 
solid  states.  Furthermore,  one  assumes  an  equal 
value  of  p°s  for  both  states.  This  implies  that  the 
radius  of  each  particle  remains  unchanged  along 
the  trajectory. 


Equilibrium  crystallization.  Assuming  an  equi¬ 
librium  liquid-to-solid  transition  process  readily  con¬ 
nects  the  specific  internal  energy,  es,  to  the  temper¬ 
ature  of  particles,  Ts,  by  the  following  relation: 


®s^s»  ^  es  <  ^s"^"m 
"^nr  ^  ^s"^"m  <  es  <  ^s"^m  + 

(es  —  hs)/Cs>  ^  ^s"*"m  +  ^*s  <  es 


Influence  on  gas.  Eqs.  (8)  and  (9)  describe  the 
changes  in  momentum  and  internal  energy  of  parti¬ 
cles  immersed  in  gaseous  flow.  The  particles’ 
effect  on  the  gas  is  taken  into  account  by  terms 
Fpart,  which  have  the  following  form: 

FPuar‘  =  2>scf(us-u) 
s 

FPvar*  =  XpsCf(vs-v)  <12) 

s 

Fhart  =  £ps{Cf[us(Us-U)  +  Vs(Vs-v)] 

+  Cq(Ts-T)  } 


where  ps  is  the  mass  density  of  the  particle  cloud. 
To  define  ps,  one  invokes  the  continuity  equation  of 
a  particle  cloud 

dPsusSs  _  0 
dt 

where  Ss  is  the  area  of  the  particle  flow-tube. 


given  previously  in  Ref.  10. 
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Numerical  Techniques 


Calculation  of  Particle  Parameters 


The  system  of  equations  describing  the  gas  and 
particle  flow  properties  is  solved  by  using  axial 
marching  numerical  methods;  that  is,  passing  from 
one  cross  section  x  =  xk  s  const  to  another,  step  by 
step  as  K  is  monitonically  increased  (k  =  0,1,2  ...). 

Passing  from  cross  section  x  =  xk  to  x  =  xk+1 
(one  step)  involves  several  iterations  which  included 
split  computations  of  gas  and  particle  parameters. 
Within  one  iteration  the  gas  parameters  are  com¬ 
puted  with  the  particle  interaction  terms,  Fparl,  taken 
from  the  previous  iteration.  Then  particle  parame¬ 
ters  are  calculated  using  the  previously  obtained 
gas  parameters.  After  that,  particle  interaction 
terms  are  recomputed  to  be  used  in  the  next  itera¬ 
tion.  In  practice,  only  two  iterations  of  this  type  are 
needed  to  define  all  gas  and  particle  properties. 

Calculation  of  Gas  Parameters 


This  section  describes  a  numerical  technique 
for  compuation  of  Eqs.  (8)-(11)  within  one  iteration 
between  the  first  (x  =  xk)  and  the  second  (x  =  xk+1) 
cross  sections. 

Let  us  consider  several  trajectories  of  the  s,h 
particle  fraction.  These  trajectories  have  numbers 
js  =  1,2 . Js.  At  the  first  section,  the  trajec¬ 

tory  coordinates,  y.k  ,  and  the  particle  parameters 
(us,  vs,  es,  r*,  fajjS  at  these  points  are  consid¬ 
ered  to  be  known.  The  gas  parameter  distributions 
are  known  at  both  x  =  xk  and  x  =  xk+1  sections. 


At  the  second  section  the  trajectory  coordinates, 
yks+ 1 ,  and  the  particle  velocity  components, 
(us,  vs)ks+ 1 ,  are  calculated  by  the  following  formu¬ 
lae: 


(vs)js  +  (vs)is+1A 

(us)k+(us)ks+1 


To  compute  the  gas  parameters  within  one  iter¬ 
ation,  the  numerical  technique  described  in  detail  in 
Ref.  1 1  is  applied.  This  technique  invokes  the  same 
assumptions  as  the  SCIPVIS  technique12  and 
allows  space-marching  calculation  of  the  subsonic 
regions  (e.g.,  subsonic  external  flow  and  flow 
regions  behind  the  Mach  disc).  An  original  second- 
order  Godunov  method  fitted  to  the  coupled  calcu¬ 
lation  of  sub-  and  supersonic  regions  is  used  to  inte¬ 
grate  gasdynamic  and  kinetic  equations  without 
using  a  splitting  technique.  An  explicit-implicit 
approximation  of  source  terms  (<5n)  on  the  right- 
hand  side  of  kinetic  Eqs.  (5)  provides  a  stable  com¬ 
putation  of  fast-relaxing  flows. 


_  (us))k8  +  0.5AtCf[Ujks  +  ufs+1(1  +  AtCQ] 
(Us)is  "  1  +AtCf(1  +0.5AtCf) 

(vs)jks  +  0.5AtCf[Vjks  +  vks+ 1  (1  +  AtC,)] 
(Vs)js  "  1 +AtCf(1 +0.5AtCf) 

where  C,  =  0.5[(Cf)k  +  (Cf)ks+1  ] ,  Ax  =  xk  +1  - xk, 

At  =  (2Ax)/[(us)k+(us)Js+1].  The  particle 
parameters  at  the  second  section,  x  =  xk+1 ,  appear¬ 
ing  in  the  right-hand  side  of  these  equations  are 
taken  from  the  previous  iteration.  The  gas  parame¬ 
ters  ^Ujk ,  uks+ 1 ,  ...j  at  both  sections  are  defined  at 
the  particle  trajectory  points. 


According  to  the  adopted  method,  the  distribu¬ 
tions  of  gas  parameters  (u,  v,  p,  P,  h0,  qn)  at  any 
cross  section  x  =  xk  are  described  by  piecewice-lin- 
ear  functions  with  discontinuities  at  the  grid  points 

y?0 =  o,  i,  **•> J;  yS < y? <  **•  < yj) *  These  9rid 

points  are  used  for  the  gas  parameter  computa¬ 
tions  only. 


In  the  part  of  the  trajectory  where  crystallization 
does  not  hold  (either  r*  =  1  and  Ts  >  Tc  or  r*  =  0) 
the  particle  temperature  is  calculated  by 
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where  Cq  =  0.5 


<°q>js  + 


In  the  part  of  the  trajectory  where  crystallization 
occurs,  Eqs.  (9)  and  (10)  are  integrated  using 
small  step  size,  At'  =  AkAt  (0  <  Ak  <  1),  succes¬ 
sively  passing  from  the  point  k'  =  k  to  the  point  k'  = 
k  +  1 ,  using  an  explicit  approach: 


,  vk  +Ak  .  .k  a  q-k  .-p  ,k 

(es)js  =  (es)js+Cq[Tjs_(Ts)jsJAt 


($*“  =  (r»)il-0T'»-(T=,!3,'84t' 

where  the  integration  step,  Af,  satisfies  the 
condition 


<T.)!i*Ak-<Tjlj'.lso-25 


Tm-<T8);; 


°4;)rksi 

In  the  presence  of  y  phase  (r*<  1  and  fa  <  1), 
Eq.  (11)  is  integrated  explicitly 


(f«)is  =  min  1,(fa)js  +  Ataaexp  -b0/(Ts)js 


Calculation  of  Gas-Particle  Interaction  Terms 

The  following  describes  an  approximation  of 
gas  particle  interaction  terms,  Fpart. 


Mass  expansion  of  the  s1h  particle  fraction 
between  the  neighboring  trajectories  is  expressed 

^  |yk 

lyjs+1 

js+1/2  - 1  J  <2*y)  Psusdv 


'JS 


The  values  of  ^js+1/2  are  specified  at  the  initial 
section,  x  =  x°,  in  compliance  with  the  problem 
conditions.  By  Eq.  (2.13),  these  values  are  not 
changed  in  the  marching  computation  progression; 
that  is,  they  are  the  same  for  any  section  x  =  xk. 

After  calculation  of  the  particle  parameters  at 
the  second  section,  x  =  xk+1,  the  following  parame¬ 
ters  describing  gas  particle  interaction  are  defined 


k  +  1/a  _  rCf(us  -  U)1 k  + 1/2  _  (us)js  -  (us)js+ 1 


U.-U 

(f  ,k+i/2=rw 

^VjS  -[  u 


Ax 


Cf(vs-v)ik  +  1/2  (vs)jks-(vs)ks+1 


JjS 


Ax 


.k.„8_fC,[U5(Us-U)tV,(V,-V)] 
js  |  Us 


Cq(Ts-T)lk  +  ,/2 _ ,,k  +  1/2 


=  (Usfu  +  Vsfv)js+ 


JS 


(es)js-(es)js 

Ax 


k+1 


where(us)ks+1^2  =0.5[(us)k  +  (us)ks+1],  (vs)k  + 1/2 
=  0.5[(vs).ks+(vs)ksk1]. 

So,  we  need  to  compute  the  gas-particle  inter¬ 
action  terms,  (Fpart)k+  \\ ,  at  the  centers  of  the 
elementary  cells  used  for  the  gas  parameter  com¬ 
putations.  Thus,  let  us  consider  the  intermediate 
section,  x  =  xk+1/2  =  0.5(xk  +  xk  +  1).  In  this  sec¬ 
tion  there  are  grid  points  yk  + 1/2  =  0.5^yk  +  y^  +  1^ 
which  define  the  intervals  (yk  +  1/2,  y’f+Y2)-  We 
will  treat  (Fpart)k+j/2  as  the  interaction  terms 
averaged  over  these  intervals. 

We  know  that  the  particle  trajectories 
vk+1/2  =  0.5fyk  +  yk+1*2l  at  the  intermediate 
section  do  not  coincide  witn  y  k  + 1/2.  Furthermore, 
the  gas-particle  interaction  terms  averaged  over 

intervals  (V[s+1/2,  yJV+i2)  can  be  aPProximated 
by  (hereafter  the  upper  index,  k  +1/2,  is  omitted) 

/pparf\  _  *15+  1/2^u)js  +  ^u^js  +  1^ 

l  u  /js  +  1/2  2[7t(yjs  +  yjs+1)]v|yjs+i  -yjs| 

/pparA  _  yjs  + 1  /2t(fy)jS  +  ^js  + 1  ] 

'  v  ^js  +  1/2  2[jt(y]s  +  yjs  +  1)]v|yjs  +  i  -yjs| 
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/'pparfj  _  Tjs  +  1/2[(fh)Js  +  ^|h^js+j2 

l  h  Jjs  +  i/2  2[jc(yjs  +  yjs  +  1)]v|yjs  +  1-yjs| 

By  this  means,  the  problem  is  how  to  compute 
(Fpart)j+i/2.  knowing  (Fpart)jS+1/2  for  each  particle 
fraction.  The  trajectories  of  two  particles  in  the 
same  size  group  may  intersect,  with  the  only 
change  being  the  reversal  of  index  of  raidal  loca¬ 
tion  for  the  two  particle. 

A  convenient  way  to  solve  this  problem  is  to 
compute  the  value  of  (Fpart)j+1/2(yj+i  -  Yj)  as  the 
sum  of  all  values  of  (Fpart)jS+i/2lyjS+i  -  Yjsl>  'n  ful1  or 
in  part,  for  which  the  interval  (yjS+i.  YjS)  is  totally 
incorporated  into  the  interval  (yj+i,  yj),  (the  value  in 
full)  or  intersects  with  it  (the  value  in  part).  If  the  lat¬ 
ter  is  the  case,  part  of  this  value  is  defined  by  the 
expression: 

(Ymax~Ymin)  x  (f')max  +  (t-)min 

|yjs+i -yjs|  (f')js+i +  (t-)js 

where  f.  =  fu,  fv  or  fh, 

Ymax  =  min[yj+i.  max  (yjs,  yjs+1)], 

ymin  =  max[yj,  min(yjs,  yjs+1)] 

(f-)max  and  (f-)min  are  calculated  at  the  points 
ymax  and  Ymin  bY  linear  interpolation  between  (f.)jS 
and  (f.)js+1. 

This  procedure  gives  a  correct  approximation  of 
(Fpart)j+1/2  for  all  cases  of  particle  trajectory  behav¬ 
iors. 

Program  Implementation 

The  described  mathematical  models  and 
numerical  techniques  have  been  embodied  in  the 
software  package,  hereafter  referred  to  as  Numeri¬ 
cal  Analysis  of  Real  Jets  (NARJ). 

NARJ  is  designed  for  the  numerical  simulation 
of  two-dimensional  steady  outflow  of  propellent 
combustion  products  from  chemical  propulsion 


devices.  NARJ  embraces  the  computations  of  gas/ 
particle  two-phase  flow  in 

•  Sub-,  trans-,  and  supersonic  parts  of  the  clas¬ 
sical  converging/diverging  nozzles; 

•  Jets  exhausted  from  the  nozzle  into  subsonic 
or  supersonic  external  gas  flow  or  into  vac¬ 
uum  free-stream  conditions. 

NARJ  consists  of  a  large  number  of  partially 
interconnected  programs  (written  in  Fortran ),  each 
solving  specific  problems  for  specified  combustion 
products.  NARJ  includes  physical-chemical 
approximations  for  chemical  reactions,  vibrational 
relaxation,  homogeneous  condensation,  and  lami¬ 
nar  and  turbulent  mixing,  and  liquid-solid  particle 
phase  change  kinetics.  Interfacing  of  the  modules 
is  required  to: 

•  utilize  output  data  from  one  program  as  an 
input  data  to  the  other  program; 

•  partially  utilize  the  same  program  modules 
when  forming  different  programs. 

Comparative  Analysis 

A  sample  calculation  was  performed,  applying 
the  NARJ  computer  program.  The  calculation  was 
initiated  with  one-dimensional  startline  conditions 
at  the  nozzle  exit  plane,  as  described  in  Tables  2 
and  3.  These  conditions  represent  a  generic,  solid- 
propellant  rocket  motor  (SRM)  with  aluminum  load¬ 
ing  of  approximately  20  percent.  The  fuel-rich  com¬ 
position  of  the  exhaust  gas,  provided  in  Table  3,  is 
also  typical  for  SRM  operation.  The  nozzle  exit 
radius,  denoted  as  Re,  was  specified  to  be  1  foot. 

In  order  to  isolate  and  compare  the  modeling 
methodologies  for  the  particle  crystallization  and 
gas/particle  coupling,  the  particles  and  gaseous 
mixture  were  assumed  to  be  initially  equilibrated, 
having  equal  temperature  and  velocity.  The  initial 
temperature  was  assumed  to  be  2,500  K,  well 
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Table  2.  One-Dimensional  Nozzle  Exit  Startline 
Conditions  for  NARJ  Plume  Simulation 


Static  Pressure 

0.4  atm  (5.9  psia) 

Static  Temperature* 

2,500  K  (4500°R) 

Axial  Velocity* 

2895  m/sec  (9498  ft/sec) 

Mach  number 

2.35 

Radial  Velocity* 

0.0  m/sec 

Al203  Density  of  each 
Size  Group 

9.07  x  10'6  gm/cm3 

Particle  Radii** 

Group  1 

1.5  pm 

Group  2 

3.0  pm 

Group  3 

6.0  pm 

*  The  gas  and  the  particulates  were  initially  assumed  to 
be  in  thermodynamic  equilibrium  having  equal  tempera¬ 
ture  and  velocity  at  the  nozzle  exit. 

**  All  particulates  are  assumed  to  be  spherical. 

Table  3.  Startline  Gaseous  Specie  Mole 
Fractions  for  NARJ  Plume 
Simulation 


CO 

0.2862 

C02 

0.01073 

Cl 

1.237E-3 

Cl2 

0.000 

H 

5.109E-3 

h2 

0.4139 

h2o 

7.409E-2 

HCI 

0.1319 

n2 

7.670E-2 

0 

1.628E-7 

OH 

6.408E-5 

02 

0.000 

above  the  assumed  NARJ  AI2O3  particle  solidifica¬ 
tion  temperature  of  2,289  K.  Therefore,  all  the  par¬ 
ticulates  were  initially  in  the  liquid  phase.  The  NARJ 
model  includes  particulate  supercooling  approxi¬ 
mation.  The  supercooling  temperature  is  specified 
as  approximately  1,900  K.  All  particle  groups  are 
expected  to  experience  the  liquid-to-solid  phase 
change  process  as  the  plume  exhaust  expands  and 
cools.  Three  particle  size  groups  were  included  in 
these  computations.  The  radii  of  the  particle  size 
groups  are  typical  for  SRM  flows  and  are  given  in 
Table  2.  The  mole  fraction  composition  of  the 
exhaust  gas  mixture  is  shown  in  Table  3.  The  mass 
percentage  of  Al203  particulates  was  approxi¬ 


mately  44  percent  of  the  total  mass,  equally  divided 
among  the  three  particle  size  groups. 

The  free-stream  ambient  flow  conditions  were 
specified  as  0.1  atm  (1.46  psia),  providing  an 
underexpanded  nozzle  operating  condition  with 
Pexit/P~=  4- The  free-stream  velocity  was  specified 
as  1,100  m/sec  (3,609  ft/sec).  The  static  tempera¬ 
ture  was  217  K  (390.6°R),  corresponding  to  a  flight 
Mach  number  of  3.7.  The  composition  of  the  ambi¬ 
ent  airstream  was  specified  by  mole  fractions  as 
[N2]  =  0.7897,  [02]  =  0.2100,  and  [C02]  =  3.0  x  10~4. 

In  this  analysis,  two  NARJ  solutions  were  com¬ 
puted.  The  first  solution  used  the  default  NARJ 
chemical  kinetic  rate  data,  based  largely  on  Rus¬ 
sian  literature  results.  For  the  second  solution, 
results  were  obtained  using  a  modified  version  of 
NARJ,  which  contained  a  set  of  chemical  reactions 
and  rates  taken  from  Baulch.13  The  intent  of  the 
second  solution  was  to  examine  the  impact  of  the 
rates  of  the  Russian  literature  as  compared  to 
those  of  the  western  literature. 


The  calculated  gas  static  temperature  contours 
obtained  from  the  default  NARJ  model  are  shown 
in  Fig.  1.  The  computational  domain  of  Fig.  1 
extends  axially  from  the  nozzle  exit  to  20  Re  down¬ 
stream.  The  atmospheric  gas  entrainment  and  mix¬ 
ing  with  the  plume  gases  contributes  to  the  shear 
layer  development,  which  is  evident  near  the  outer 
periphery  of  the  contour  map.  The  shear  layer  is 
delineated  from  the  inviscid  core  flow  by  the  tem- 


0  20  Re 

Fig.  1 .  Spatial  map  of  default  NARJ  solution  for  the 


gas  temperature. 
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perature  gradient  showing  a  region  of  elevated 
temperatures  in  the  outer  radial  expansion  region. 
Combustion  and  shear  heating  contribute  to  the 
temperature  gradient  in  the  shear  layer.  The  one¬ 
dimensional,  inviscid  core  flow  is  evident  near  the 
nozzle  exit,  and  persists,  undisturbed,  along  the 
centerline  for  approximately  1  nozzle  radius. 

Figure  2  is  the  static  pressure  contour  map  for 
the  same  nearfield  segment  of  the  NARJ  computa¬ 
tion.  The  rapid  expansion  of  the  underexpanded 
plume  exhaust  near  the  nozzle  lip  is  evident.  An 
inviscid  wave  structure  emanating  in  the  near-field 
region  is  evident  and  can  be  tracked  as  it  propa¬ 
gates  and  reflects  off  the  plume  centerline  and  the 
plume  outer  boundary.  This  wave  propagation  is 
also  observed  in  the  gas  temperature  map  of  Fig.  1 
as  elevated  temperatures  along  the  centerline,  10 
Re  downstream. 


0  20  Re 

Fig.  2.  Spatial  map  of  default  NARJ  solution  for 
the  logarithm  of  the  gas  pressure  in 
Pascals. 

Fig.  3  is  the  calculated  contour  map  of  H20  mole 
fraction  showing  the  kinetic  rate  determination  of 
H20  concentration  using  the  NARJ  default  chemis¬ 
try.  The  factor  of  2  increase  in  the  H20  mole  fraction 
evident  in  comparing  the  initial  H20  mole  fraction 
level  to  the  shear  layer  H20  mole  fraction  value  indi¬ 
cates  significant  combustion  resulting  in  formation 
of  H20  in  the  shear  layer.  The  formation  of  H20  in 
the  shear  layer  results  from  ambient  02  entrain¬ 
ment,  mixing,  and  combustion  with  the  excess  H2  in 
the  plume  exhaust  gas  mixture.  It  is  interesting  to 
note  the  very  gradual  initial  development  of  the 


shear  layer  and  the  accompanying  combustion 
zone.  Initially,  there  is  very  little  production  of  H20 
near  the  nozzle  lip,  where  the  ambient  air  is  not  suf¬ 
ficiently  entrained  and  the  mixing  is  slight. 


o  20  Re 

Fig.  3.  Spatial  map  of  default  NARJ  solution  for 


the  H20  mole  fraction. 

Figure  4  is  the  centerline  axial  profile  of  H2, 
H20,  and  N2  mole  fractions,  showing  the  creation 
of  H20  and  the  corresponding  H2  depletion  and  the 
mixing  of  the  atmospheric  N2  with  the  plume 
gases.  The  centerline  static  temperature  of  the  gas 
is  also  included  on  Fig.  4.  The  afterburning  shear 
layer  appears  to  reach  the  plume  centerline  posi¬ 
tion  at  roughly  50  Re,  as  indicated  by  the  presence 
of  the  N2  from  the  ambient  air,  and  extends  consid¬ 
erably  into  the  far  field,  to  roughly  300  Re.  The 
increase  in  the  temperature,  followed  by  the  grad¬ 
ual  decrease  in  this  region,  is  also  indicative  of 
shear  layer  combustion.  Downstream  of  the  loca¬ 
tion  of  air  entrainment  on  the  centerline,  the  core 


Fig.  4.  Axial  profiles  of  the  default  NARJ  jet  center- 
line  properties. 
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flow  field  is  dominated  by  mixing  and  combustion 
until  the  excess  H2  is  depleted  near  200  Re. 

As  suggested  in  the  previous  description  of  the 
NARJ  model,  simulation  of  the  coupled  Al203  par¬ 
ticulate  properties,  including  drag,  trajectory,  and 
phase  change  thermodynamics  constitutes  an 
essential  portion  of  the  computation  methodology. 
Near-field  particle  temperature  contours  for  the  1 .5- 
pm  radius  particle  group  is  shown  in  Fig.  5.  This 
contour  map  focuses  on  the  near  field  region  of  the 
computational  domain,  extending  axially  to  20  Re. 
The  small  particulates  experience  the  smallest  drag 
force  per  particle  and  are  nearly  in  equilibrium  with 
the  gas.  As  such,  small  particulates  display  a  fairly 
rapid  cooling  trend  similar  to  the  gas  phase  temper¬ 
ature  trend  (shown  previously  in  Fig.  1).  The  small 
particles  track  the  gas  flow  in  the  shear  layer  and 
are  heated  by  the  combustion  taking  place  in  this 
region.  The  particle  temperatures  also  increase  in 
the  downstream  region  near  the  centerline.  This 
sudden  temperature  rise  near  the  centerline  is  due 
to  the  phase  change  “blik”  phenomenon  that  results 
when  the  supercooled  liquid  particulates  rapidly 
transition  to  the  metastable  gamma  solid  phase, 
accompanied  by  a  rapid  rise  in  temperature.  This 
phenomenon  is  also  evident  in  Fig.  6,  a  contour 
map  of  the  solid  metastable  gamma  crystalline 
phase  fraction  for  the  1.5-|im  Al203  particulates. 

The  liquid  phase  is  prevalent  in  the  near-field 
region  of  the  jet,  as  expected  since  the  particles 
were  initiated  in  the  liquid  phase.  The  first  occur¬ 
rence  of  the  metastable  gamma  solid  phase  is 
downstream  near  the  centerline,  where  the  super- 


0  20  Re 

Fig  5.  Spatial  map  of  default  NARJ  solution  for  the 
temperature  of  the  1 .5-p.m  particle  group. 


0  20  Re 

Fig.  6.  Spatial  map  of  default  NARJ  solution  for 

the  gamma  phase  fraction  of  the  1 .5-pm 

particle  group. 

cooling  temperature  “blik”  phenomenon  occurs.  As 
soon  as  the  particle  supercooling  temperature  is 
realized,  the  transition  process  from  liquid  to  meta¬ 
stable  gamma  crystalline  structure  begins,  coinci¬ 
dent  with  the  kinetically  controlled  transition  from 
the  metastable  gamma  to  the  stable  alpha  crystal¬ 
line  structure.  The  spatial  distribution  of  the  alpha 
phase  fraction  for  this  particle  group  is  shown  in  Fig. 
7.  As  expected,  the  alpha  phase  formation  occurs 
within  the  same  spatial  region  as  the  gamma  phase. 
The  alpha  phase  production  rate  is  extremely  tem¬ 
perature  sensitive,  decreasing  by  nearly  6  orders  of 
magnitude  as  the  temperature  declines  from  the 
solidification  temperature  of  2,289  K  to  1,500  K. 
Therefore,  the  majority  of  the  formation  of  alpha 
phase  occurs  in  the  spatial  region  immediately  fol¬ 
lowing  the  particle  “blik,”  where  the  particles  are  ini¬ 
tially  solidified  at  relatively  high  temperatures.  Over¬ 
all,  the  crystalline  phase  transition  process  does  not 
appear  to  influence  the  gas  dynamic  behavior  of 
any  portion  of  the  flow.  However,  differences  in  the 
solid  phase  crystalline  structure  could  influence  the 
radiative  transfer  characteristics  of  the  jet  if  the 
gamma  and  alpha  crystalline  structures  have  signif¬ 
icantly  distinct  optical  properties. 


2.0E-03 

1.6E-03 

1.2E-03 

8.0E-04 

4.0E-04 

0.0E+00 


0  20  Re 

Fig.  7.  Spatial  map  of  default  NARJ  solution  for  the 
alpha  phase  fraction  of  the  1.5-pm  particle 
group. 
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The  two  larger  particle  groups  experience 
greater  drag  force  per  particle,  resulting  in 
decreased  cooling  rates  and  velocities.  Therefore, 
the  particle  phase  change  solidification  process  is 
pushed  farther  downstream  beyond  20  Re  axial 
distance. 

The  NARJ  model  was  modified  to  replace  the 
Russian-based  chemical  reaction  rates  with  a  com¬ 
mon  set  of  chemical  equations  and  rates  based  on 
Baulch,  et  al.13  The  temperature-dependent  rate 
equations  and  the  associated  rate  coefficients  are 
tabulated  in  Table  4. 

A  comparison  of  the  rate  coefficients  in  Table  4 
with  those  of  Table  1  indicates  that  most  rates  are 
similar;  however,  some  differences  are  noted.  The 
most  obvious  is  the  exclusion  of  the  nitrogen  reac¬ 
tion  set  in  the  modified  chemical  equation  set.  Fur¬ 
thermore,  the  reactions  leading  to  the  exothermic 
formation  of  H20  and  C02  are  considerably  faster 
in  the  default  NARJ  chemical  rate  system.  The 
exception  to  this  is  the  CO  +  O  +  M  -» C02  +  M  reac¬ 
tion,  in  which  the  default  NARJ  rate  is  slower. 

The  global  effect  of  the  different  kinetic  reaction 
data  is  noted  by  contrasting  NARJ  solutions  incor- 
Table  4.  Chemical  Reaction  Rate  Coefficients  used 


the  Modified  NARJ  Computations  K+  =  AT  n 
exp(-E/RT) _ _ _ _ 


Reaction 

A,  m6/(kmol-sec)2 

n 

E,  cal 

h  +  h  +  m  =  h2  +  m 

1.09E+12 

1 

0 

H  +  0  +  M  =  OH  +  M 

3.63E+12 

1 

0 

0  +  0  +  M=  02  +  M 

1 .09E+08 

0 

-1,800 

H  +  OH  +  M  =  HaO  +  M 

1.09E+08 

2 

0 

CO  +  0  +  M  =  C02  +  M 

1 .09E+08 

0 

4,360 

Cl  +  Cl  +  M  =  Cl2  +  M 

2.36E+09 

0 

-1,800 

H  +  Cl  +  M  =  HCI  +  M 

3.63E+08 

2 

0 

Reaction 

A,  m3/kmol-sec 

n 

E,  cal 

OH  +  0  =  H  +  02 

1.81E+10 

0 

960 

OH  +  H  =  H2  +  0 

8.43E+06 

-1 

7,000 

OH  +  OH  =  H20  +  0 

6.02E+09 

0 

1,100 

H2  +  OH  =  H20  +  H 

2.11E+10 

0 

5,180 

oh+co  =  co2  +  h 

1 .69E+04 

-1.3 

-660 

H20  +  Cl  =  HCI  +  OH 

8.43E+09 

0 

4,260 

HCI  +  0  =  OH  +  Cl 

3.01E+10 

0 

17,600 

H  +  Cl2=  HCI  +  Cl 

3.61  E+09 

0 

6,000 

H  +  Cl2  =  HCI  +  Cl 

9.03E+10 

0 

1,200 

porating  the  default  chemistry  with  the  version  of 
NARJ  incorporating  the  modified  rates.  Figure  8a 
contrasts  centerline  axial  profiles  of  the  gas  static 
temperature  from  NARJ  solutions  obtained  with 
default  and  modified  rates.  The  immediate  near 
field  results  are  identical,  as  the  chemistry  is  essen¬ 
tially  frozen  until  the  mixing  layer  intersects  the 
plume  centerline  at  an  axial  distance  of  approxi¬ 
mately  50  Re,  denoting  the  beginning  of  the  center- 
line  afterburning  region.  The  location  where  the 
mixing  layer  reaches  the  centerline  is  clearly  evi¬ 
dent  in  Fig.  8b,  which  compares  the  axial  profiles  of 
several  chemically  active  species  obtained  from  the 
default  and  the  modified  chemistry  models.  Based 
upon  the  depletion  of  H2  and  the  generation  of  H20, 
the  afterburning  region  extends  to  approximately 


a.  Gas  temperature 

— Def  Rates:  H20  — * —  Def  Rates:  OH 

_ Mod  Rates:  H20 Mod  Rates:  OH 

A  —  Def  Rates:  C02  — ¥ —  Def  Rates:  H2 
—A—  Mod  Rates:  C02  Mod  Rates:  H2 


b.  Gas  composition 

Fig.  8.  Axial  profiles  on  the  jet  centerline  using  the 
default  NARJ  rates  and  the  modified  rates. 


12 

American  Institute  of  Aeronautics  and  Astronautics 


Fig.  9.  Predicted  NARJ  pressure  map.  The  values  are  the  logarithm  of  the  pressure  in  atm. 


300  Re  in  both  simulations.  However,  the  NARJ 
default  chemistry  generates  consistently  higher  gas 
temperatures.  Downstream  of  the  afterburning 
region,  diffusive  cooling  dominates  the  flow,  and 
the  two  solutions  approach  the  same  temperature 
values.  Thus,  while  the  modified  rates  do  signifi¬ 
cantly  alter  the  gas  temperature  in  the  afterburning 
region  and  slightly  modify  the  species  concentra¬ 
tions,  the  different  chemical  equations  and  rates  do 
not  significantly  alter  the  spatial  location  of  the  tem¬ 
perature  maximum,  nor  the  global  structure  and 
composition  of  the  non-afterburning  jet  regions.  For 
the  sake  of  this  comparative  analysis,  all  subse¬ 
quent  NARJ  solutions  were  obtained  incorporating 
the  modified  reactions  and  rate  coefficients. 

In  order  to  further  assess  the  NARJ  model,  a 
computed  near-field  pressure  map  is  shown  in  Fig. 
9.  The  axial  extent  of  Fig.  9  is  46  Re.  The  nozzle 
exit  radius  is  1  foot.  The  intensity  values  are  dis¬ 
played  as  the  logarithm  of  the  pressure  in  atmo¬ 
spheres.  The  NARJ  solution  indicates  a  wave  prop¬ 
agation  structure  in  the  near  field.  Following  the  ini¬ 
tial  expansion  region,  the  pressure  is  nearly  equal 
to  the  ambient  pressure  in  the  remaining  computa¬ 
tional  regions.  The  pressure  downstream  of  5  Re,  is 


Fig.  10.  NARJ  axial  profiles  of  the  gas  pressure  on 
the  jet  centerline. 

nearly  equilibrated  with  the  ambient  pressure  (Fig. 
10)  which  displays  the  centerline  axial  gas  static 
pressure  profile.  The  small  spike  near  2  Re  is  evi¬ 
dence  of  a  pressure  wave  caused  by  the  initiation 
of  burning  in  the  shear  layer  immediately  down¬ 
stream  of  the  nozzle  exit.  Following  the  precipitous 
pressure  drop,  the  NARJ  centerline  result  oscillates 
about  the  ambient  pressure. 

Figure  11  shows  the  NARJ  gas  temperature 
contour.  The  axial  extent  of  the  temperature  con¬ 
tour  map  is  46  Re.  Note  that  the  NARJ  radial  tem¬ 
perature  profiles  have  a  local  minimum  on  the  cen- 


Fig.  1 1 .  Predicted  NARJ  temperature  map. 


2.60E+03 

2.42E+03 
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terline  with  a  local  maximum  at  the  shear  layer  for 
axial  distances  less  than  50  Re. 

The  centerline  temperature  profiles  of  gas  and 
particulates  are  shown  in  Fig.  12.  The  relatively  low 
NARJ  gas  temperature  allows  the  solid  phase  tran¬ 
sition  to  occur  in  the  smallest  particle  groups  near 
the  nozzle  exit  plane,  as  the  temperature  drops 
below  1 ,900  K  in  this  region.  Thus,  in  this  simula¬ 
tion,  both  the  1 .5-|xm  radius  and  the  3.0-pm  radius 
particle  groups  undergo  liquid  to  solid  transitions 
prior  to  the  occurrence  of  the  afterburning  region. 
The  current  NARJ  model  does  not  incorporate  a 
mechanism  allowing  the  solid  particulates  to  re¬ 
melt  to  the  liquid  phase,  a  phenomenon  which 
could  potentially  occur  in  downstream  combustion 
regions.  The  absence  of  the  second  phase  change 
process  could  artificially  enhance  the  proportion  of 
the  alpha  solid  phase  in  the  downstream  region. 


Fig.  12.  NARJ  axial  profiles  of  the  gas  and  parti¬ 
cle  temperatures  on  the  jet  centerline. 

Figure  13  shows  the  centerline  temperature  pro¬ 
files  of  the  NARJ  particles  and  gas  for  the  complete 
solution  extending  to  500  Re,  while  Figure  14  is  the 
centerline  axial  profiles  of  the  mole  fractions  of  N2 
and  OH.  The  axial  location  for  the  appearance  of 
N2  is  in  keeping  with  the  axial  location  of  the  center- 
line  temperature  rise,  indicating  the  first  occurrence 
of  the  shear  layer  mixing  at  the  centerline.  This 
temperature  rise  delays  the  phase  transition  of  the 
largest  particle  group  to  downstream  of  the  after¬ 
burning  region. 


particle  temperatures  on  the  jet  center- 
line,  using  the  modified  NARJ  rates. 


Fig.  14.  NARJ  axial  profiles  of  the  OH  and  N2  mole 
fractions  on  the  jet  centerline,  using  the 
modified  NARJ  rates. 

Conclusions 

This  work  describes  and  demonstrates  a  com¬ 
prehensive  computational  methodology  for  simulat¬ 
ing  the  jet  flow-field  properties  of  multiphase,  chem¬ 
ically  reacting,  aluminum-loaded  propulsion-gener¬ 
ated  exhaust.  The  NARJ  model  incorporates  a  wide 
range  of  physical  phenomena,  including  vibrational 
relaxation  and  vapor  condensation,  which  were  not 
explicitly  discussed  in  this  work.  Physical  approxi¬ 
mations  incorporated  into  the  model  include  cou¬ 
pled  temperature-dependent  chemical  kinetic 
rates,  two-phase  gas/particle  drag,  and  turbulent 
mixing  and  particle  phase  change  thermodynamics 
including  supercooling  phenomena.  The  NARJ 
methodology  solves  the  inviscid  and  viscous  flow 
regions  using  a  unified  approach.  It  also  includes  a 
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two-step  kinetic  AI2O3  particle  phase  change  model 
in  which  the  phase  change  process  includes  an  ini¬ 
tial  transition  from  liquid  to  solid  metastable  gamma 
crystalline  structure  immediately  followed  by  a 
kinetic  controlled  rate  transition  to  the  stable  alpha 
solid  phase  crystalline  form.  The  transition  from 
gamma  to  alpha  is  temperature  dependent,  occurs 
rapidly  at  elevated  temperatures,  and  is  extremely 
slow  at  moderate  temperate  ranges.  Therefore,  the 
majority  of  the  solid  AI2O3  in  jet  flows  is  predicted  by 
NARJ  to  be  in  the  metastable  gamma  phase. 
These  studies  conclude  that  the  crystalline  struc¬ 
ture  formation  does  not  significantly  influence  the 
gas  dynamic  properties;  however,  it  could  be  an 
issue  for  radiative  transfer  simulations  if  the  gamma 
and  alpha  phases  have  substantially  different  opti¬ 
cal  properties. 
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